/* ========================================================================
 * sim_mc.c
 * The simulator implementation of Markov chain models
 * 
 * Andr� Pedro (c) apedro@gmx.com
 * 
 * Discrete event system simulator based on a Markov chain model
 *   Simulates a generelized semi-Markov process
 *
 * ======================================================================*/
/* $Revision: 1 $ */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#include "sim_mc.h"
#include "sim_interface.h"
#include "fibheap.h"
#include "mt19937.h"

#ifdef MATLAB_MEX_FILE   
#include "mex.h"
#define debug_error mexPrintf
#else
#define debug_error printf
#endif

/* Generate a exponential random value from the inverse function method */
double expd_sample(double lambda) {
    double U;
    
    /* Change random generation to more fast version implementation*/
    /*U = (double)rand() / ((double)RAND_MAX + 1);*/
    U = genrand_real3();
    
    return(-lambda*log(1-U));
}

/* Generate a weibull random value from the inverse function method */
double wbd_sample(double lambda, double k) {
    double U;
    
    if (!(k>0))
        return 0.0;
    
    /* Change random generation to more fast version implementation*/
    /*U = (double)rand() / ((double)RAND_MAX + 1);*/
    U = genrand_real3();
    
    return (lambda*pow(-log(1-U), 1.0/k));
}

double logn_sample(double mu, double sigma) {
    double U1, U2, v1, v2, w, Y;
    
    while(1) {
        U1 = genrand_real3();
        U2 = genrand_real3();
        
        v1 = 2*U1-1;
        v2 = 2*U2-1;
        w = v1*v1 + v2*v2;
        if(w < 1.0) {
            Y = pow((-2*log(w))/w, 1/2);
            return exp(mu + sigma * v1 * Y);
        }
    }
}

void sim_mc(const int nsim, double seed, double *S, const int n_S,
            double *E, const int n_E, double *ES,
            double *p, dparam *dp, double *C, int *path) {
    
    int    i, j, e;
    int    n, s, eid_val;
    unsigned long nseed, active_events;
    double sample, old_winner_sample, eclk_val;
    /* The fibonaci heap*/
    fibHeap* fh;
    fibNode **schevents;
    fibNode *tmp_fibnode;
    fh = make_fibheap(); /* Heap initialization */
    
    /* Change random generation to more fast version implementation*/
    /* Initialization of the random generator */
    nseed = (unsigned long)(time(NULL)*seed);
    init_genrand(nseed);
    /*DEBUG mexPrintf("Seed: %i\n", nseed);*/
    /*srand(time(NULL));*/
    
    /* Initialization */
    n=1;
    s=1;
    eclk_val=0.0;
    eid_val=0;
    
    /* Initialize schevents array */
    /* Allocate references to every node of the heap that is stored
     * in a schevents vector */
    schevents = malloc(sizeof(fibNode *)*n_E);
    
    for (n=1; n <= nsim; n++) {
        /* Set all clocks to active events in state s,
         * and the inactive events to a high clock value */
        for (i=1, active_events=0; i <= n_E; i++) {
            /*DEBUG debug_error("ES index: %d with value %d\n",(s-1)+((n_S)*(i-1)), ES[(s-1)+((n_S)*(i-1))]);*/
            /* Active clock's e[i] to state s*/
            if (ES[(s-1)+((n_S)*(i-1))]) {
                active_events++;
                /* Select sample from adequate distribution function (Exponential distribution; Weibull distribution)*/
                if (strcmp(dp[i-1].nameid, "exp") == 0) {
                    sample = expd_sample(dp[i-1].param->lastNode->data);
                    /*DEBUG mexPrintf("Expd(%f): %f\n", dp[i-1].param->lastNode->data, sample);*/
                }
                else if (strcmp(dp[i-1].nameid, "wbl") == 0) {
                    sample = wbd_sample(dp[i-1].param->lastNode->data, dp[i-1].param->lastNode->prev->data);
                    /*DEBUG mexPrintf("Wbd(%f,%f): %f\n", dp[i-1].param->lastNode->data, dp[i-1].param->lastNode->prev->data, sample);*/
                }
                else if (strcmp(dp[i-1].nameid, "logn") == 0) {
                    sample = logn_sample(dp[i-1].param->lastNode->data, dp[i-1].param->lastNode->prev->data);
                    /*DEBUG mexPrintf("Logn(%f,%f): %f\n", dp[i-1].param->lastNode->data, dp[i-1].param->lastNode->prev->data, sample);*/
                } else
                    mexErrMsgTxt("PDF not match.");
                
                if (n==1) {
                    /* set initialized clock to the event scheduller (used to select the triggered event) */
                    schevents[i-1] = make_fibnode(sample, i);
                    fibheap_insert(fh, schevents[i-1]);
                }
                else if (n>1) {
                    /* Verify if is an old event or a new active event*/
                    if (schevents[i-1]->key == (double) INT_MAX) {
                        /* New event ocurred */
                        fibheap_decreasekey(fh, schevents[i-1], sample);
                    }
                    else {
                        /* Old event ocurred */
                        /* sample value = key - old_winner_sample*/
                        sample = schevents[i-1]->key - old_winner_sample;
                        fibheap_decreasekey(fh, schevents[i-1], sample);
                    }
                }

                /* log clock values in variable C whith the calculated clock samples */
                C[(n-1)+((nsim)*(i-1))] = sample; /* set this to exponential distribution */
                /*mexPrintf("%d %f %f-%f\n", (n-1)+((nsim)*(i-1)), sample);*/
            }
            /* Inactive clock's e[i] to state s*/
            else if (!ES[(s-1)+((n_S)*(i-1))]) {
                /* Set all inactive clock events e[i] to high value*/
                if (n==1) {
                    /* set initialized clock infinitelly to the event scheduller (used to set an inactive events) */
                    schevents[i-1] = make_fibnode((double) INT_MAX, i);
                    fibheap_insert(fh, schevents[i-1]);
                }
                else if (n>1) {
                    /* verifies if a clock is a new clock or an old clock*/
                    if (schevents[i-1]->key != (double) INT_MAX) {
                        /* remove node*/
                        fibheap_decreasekey(fh, schevents[i-1], -1.0);
                        tmp_fibnode = fibheap_extractmin(fh);
                        /* insert as unitialized clock*/
                        tmp_fibnode->key = (double) INT_MAX;
                        fibheap_insert(fh, tmp_fibnode);
                    }
                }
                
                C[(n-1)+((nsim)*(i-1))] = (double) -1.0; /* Set this to max double */
            }
        }
        
        /* If all events are unavailable*/
        if(active_events == 0)
            break;
        
        /* Triguer the event with the small clock value */
        tmp_fibnode = fibheap_extractmin(fh);
        old_winner_sample = eclk_val = tmp_fibnode->key;
        eid_val = tmp_fibnode->id;
        
        /* Set winner clock to high value*/
        tmp_fibnode->key = (double) INT_MAX;
        fibheap_insert(fh, tmp_fibnode);
        
        /* Save the path (S,E) */
        path[(n-1)] = s;
        path[(n-1)+(nsim)] = eid_val;

        /* Save the hold time in C matrix */
        C[(n-1)+((nsim)*n_E)] = eclk_val;

        /* Select the next state; P(next state | actual state, event) */
        j=s-1;
        e=eid_val-1;
        for(i=0; i<n_S; i++) {
            if (p[(i+(j*n_S))+(n_S*n_S*e)] == 1)
                s=i+1;
        }
        
        /*DEBUG printf("Next state: %d\n", s);*/
        
    }
    
    /*free(clk);*/
    /*fh_deleteheap(clk_heap);*/
}
